function df=NumJacob(f,x0,varargin)

epsilon = 1e-6;
l_x0=length(x0);
f0=feval(f,x0,varargin{:});
l_f=size(f0,1);
for i=1:l_x0
    dx=[zeros(i-1,1);epsilon;zeros(l_x0-i,1)];
    df(:,i)=(feval(f,x0+dx,varargin{:})-f0)/epsilon;
end